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NUMERICAL SOLUTIONS OF THE NAVTER-STOKES EQUATIONS 
FOR THE SUPERSONIC LAMINAR FLOW OVER A 
TWO-DIMENSIONAL COMPRESSION CORNER* 

By James E. Carter 
Langley Research Center 

SUMMARY 

Numerical solutions have been obtained for the supersonic, laminar flow over a 
two-dimensional compression corner. These solutions were obtained as steady -state 
solutions to the unsteady Navier -Stokes equations using the finite -difference method of 
Brailovskaya, which has second-order accuracy in the spatial coordinates. Good agree- 
ment was obtained between the computed results and the wall pressure distributions mea- 
sured experimentally by Lewis, Kubota, and Lees for Mach numbers of 4 and 6.06, and 
respective Reynolds numbers, based on free -stream conditions and the distance from the 
leading edge to the corner, of 6.8 x 10^ and 1.5 x 10^. In those calculations, as well as in 
others, sufficient resolution was obtained to show the streamline pattern in the separation 
bubble. Upstream boundary conditions to the compression -corner flow were provided by 
numerically solving the unsteady Navier-Stokes equations for the flat -plate flow field, 
beginning at the leading edge. The compression-corner flow field was enclosed by a 
computational boundary with the unknown boundary conditions supplied by extrapolation 
from internally computed points. Numerical tests were performed to deduce that the 
magnitude of the errors introduced by the extrapolation was small. 

Calculations were made to show the effect of ramp angle and wall suction on the 
interaction flow field. The pressure distributions obtained in the present calculations, 
including a case of incipient separation, were plotted together by using the free -interaction 
scaling of Stewartson and Williams. A good correlation of the numerical results was 
found, but only fair agreement was found between this correlation and the universal pres- 
sure distribution found numerically by Stewartson and Williams. 


* The information presented herein is based in part upon a thesis entitled "Numeri- 
cal Solution of the Supersonic, Laminar Flow Over a Two-Dimensional Compression 
Corner" submitted in partial fulfillment of the requirements for the degree of Doctor of 
Philosophy in Aerospace Engineering, Virginia Polytechnic Institute and State University, 
Blacksburg, Virginia, August 1971. 



INTRODUCTION 


A problem that has interested fluid dynamicists for a number of years is the super- 
sonic flow over a flat plate followed by a ramp. The pressure rise generated by the ramp 
extends upstream along the flat plate and results in a complex interaction between the 
boundary layer and the outer inviscid stream. This interaction leads to flow separation 
for certain ranges of the Mach number, Reynolds number, and ramp angle. The purpose 
of the present investigation was to obtain numerical solutions to the finite -difference form 
of the Navier-Stokes equations for the laminar flow over a two-dimensional compression 
corner. In addition to its theoretical interest, this problem is of practical importance in 
predicting the pressure and heat loads in a wing -flap juncture on a supersonic aircraft. 
When flow separation occurs, reduced flap effectiveness results, and in the reattachment 
region the surface heating may become severe. 

A schematic diagram of the supersonic flow over a compression corner is shown 
in figure 1. The adverse pressure gradient generated by the ramp thickens the approach- 
ing boundary layer. The inner part of the boundary layer near the surface may have an 
insufficient total pressure and thus may separate from the surface since it cannot over- 
come the adverse pressure gradient. The separated boundary layer now becomes a free 
shear layer external to a steady, recirculating inner flow near the corner. The boundary 
between the shear layer and recirculation region is logically referred to as a dividing 
streamline. Farther downstream the shear layer impinges on the ramp in the reattach- 
ment region; the flow in the boundary layer then continues to accelerate until the boundary 
layer reaches a minimum thickness at the "neck." Downstream of the neck, the boundary 
layer returns to a normal state of weak interaction at the new Mach number. 

Reattachment 
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A direct method of including the upstream influence in flow fields involving large 
interaction between the viscous and inviscid flow is to treat the full Navier-Stokes equa- 
tions. Treatment of the Navier-Stokes equations avoids the uniqueness questions inherent 
in the boundary -layer interaction approach. The accuracy is probably improved since, 
for example, a solution of the Navier-Stokes equations is required in the immediate vicin- 
ity of a sharp corner as mentioned by Van Dyke (ref. 1). In addition, numerical solutions 
to the Navier-Stokes equations, although still limited by computer size and speed, serve 
as benchmark solutions for comparison with approximate but more expedient methods. 

In the present investigation, the unsteady Navier-Stokes equations are solved by 
the explicit finite -difference scheme of Brailovskaya (ref. 2). Including the unsteady 
terms in the Navier-Stokes equations results in a parabolic system of partial differential 
equations and in a well-posed, initial -value — boundary -value problem. As indicated by 
Crocco (ref. 3), the inclusion of time in the equations allows the solution to progress 
naturally from an initial guess to an asymptotic state, which is the solution to the steady 
equations. In the present investigation, it was assumed that a steady solution exists for 
the compression- corner flow field. 


BACKGROUND 
Boundary -Layer Methods 

Previous theoretical treatments of this problem have been made with the boundary - 
layer equations to describe the separation of a supersonic flow, and these methods are 
applicable to the compression-corner flow field. Three of these methods are capable of 
obtaining a solution in the reattachment region as well as in the separation region. These 
three methods are the integral methods of Lees and Reeves (ref. 4), and Nielsen, Goodwin, 
and Kuhn (ref. 5), and the finite -difference method of Reyhner and FKIgge-Lotz (ref. 6). 
Although these methods differ in detail, their success in this problem is based on several 
features which they have in common, and hence they will be discussed collectively. 
Numerical comparisons of these methods with available experimental data have been 
made by Murphy (ref. 7) and, to a lesser extent, by Hill (ref. 8). 

Standard procedures for solving the boundary -layer equations utilize forward- 
marching techniques with a known pressure distribution. Such a procedure is not appli- 
cable to the compression-corner flow field since the pressure distribution is not known 
a priori. As a result, these methods make repeated iterations in order to obtain a unique 
solution which properly accounts for the upstream influence of the corner. 

This iterative procedure is initiated at an arbitrary point upstream of the corner by 
introducing either a small increment in the pressure, as in the methods of Lees and 
Reeves and that of Nielsen et al. or by subjecting the boundary layer to a small adverse 



pressure gradient, as in the method of Reyhner and Fliigge-Lotz. The resulting pressure 
distribution is computed from the Prandtl-Meyer relation which relates the pressure to 
the local flow inclination at the boundary -layer edge. Both techniques initiate an amplify- 
ing process. As the pressure increases, the boundary -layer thickness and slope increases, 
and, in turn, causes the pressure to increase further by an amount given by the Prandtl- 
Meyer equation. This process culminates in separation of the flow which relieves the 
pressure gradient, and therefore the downstream flow approaches a constant -pressure 
region (a plateau). 

Once the plateau region is reached, the disturbance that induced the separation 
(for example, an incident shock or compression ramp) is introduced. An infinite family 
of solutions can be found by varying the strength of the disturbance (alternately, the dis- 
turbance strength can be held constant and the upstream point of interaction varied). 
Uniqueness is established by imposing downstream compatibility conditions and solving 
iteratively until these conditions are satisfied. The downstream compatibility condition 

in the methods of Nielsen et al. and Reyhner and Fliigge-Lotz is that =0 at 

dx ^2 

the same point, whereas in the method of Lees and Reeves the downstream solution is 
required to approach the flat -plate solution. 

A somewhat different approach has been taken by Stewartson and Williams (ref. 9), 
who found a universal solution to the boundary -layer equations for the region from the 
start of the interaction to a point just downstream of separation. The existence of a 
universal solution was suggested by the pressure correlation discovered experimentally 
by Chapman, Kuehn, and Larson (ref. 10) and has been verified more recently by Lewis, 
Kubota, and Lees (ref. 11). They observed that in a supersonic flow field, the initial 
pressure rise through separation to the plateau value was independent of the details of the 
disturbing mechanism (and hence is referred to as a "free interaction"), whether it be, 
for example, an incident shock, a forward -facing step, or a ramp. Stewartson and 
Williams use the method of matched asymptotic expansions to show that the length scale 
of the free -interaction region is R^xq w ^ ere ^xq is the Reynolds number based on 
free -stream conditions and the length from the leading edge to the start of the interaction. 
This scaling was the same as that found earlier by Gadd (ref. 12) by a more approximate 
method. From the expansion procedure it becomes apparent that an inner boundary layer 
of constant density and large velocity perturbation is the key feature of the free -interaction 
zone. The idea of an inner boundary layer, that is, a region where the disturbances to the 
viscous forces are comparable with the disturbances to the inertial forces, was originated 
by Lighthill (ref. 13) in studying the interaction between the boundary layer and a weak 
shock. Stewartson and Williams obtained a universal pressure distribution for the free- 
interaction region by solving the incompressible boundary -layer equations for the inner 


#* 


4 



region subjected to novel boundary condition. Only fair agreement was obtained between 
their theoretical result and an experimental wall pressure distribution obtained by 
Chapman et al. (ref. 10). 

Numerical Solution of Navier-Stokes Equations 

In the last several years the volume of literature has increased rapidly on the solu- 
tion of various problems by treating the finite -difference form of the unsteady Navier- 
Stokes equations. Recently, Cheng (ref. 14) reviewed the literature and reported most of 
the investigations which have been made. Several numerical investigations have been 
made for the flow over a rearward-facing step. Allen and Cheng (ref. 15) obtained numer- 
ical solutions for this flow field for a Reynolds number, based on the base half-height and 
the inflow conditions, of less than 1000, and a Mach number range from 2 to 4. The entire 
flow field was enclosed by a computational boundary, and the flow conditions were either 
known or had to be approximated along this boundary. Allen and Cheng used a modification 
of the explicit, time -dependent method introduced by Brailovskaya (ref. 2). Because of 
the low Reynolds number and low value of density which occurs in the recirculation region 
behind the step, they modified the Brailovskaya scheme to eliminate the dependence of the 
time step on the kinematic viscosity coefficient. This work has been extended by Ross and 
Cheng (ref. 16) to include variable viscosity and base injection. The results obtained by 
Allen and Cheng and by Ross and Cheng appear to be qualitatively correct since the char- 
acteristic features of the base flow region are evident. The quantitative results are yet 
to be verified by experiment since no experiments have been performed for the low 
Reynolds numbers used in these calculations. 

Roache and Mueller (ref. 17) reported calculations for both compressible and incom- 
pressible flow over a rearward -facing step. They used the first-order windward differ- 
ence scheme which is known to suffer from a diffusive truncation error. Unfortunately, 
it is difficult to assess the effects of this truncation error in a complicated flow -field 
calculation. Carter (ref. 18) made numerical studies of Burger's equation, which is a 
simple one -dimensional model of the Navier-Stokes equations, and showed that the wind- 
ward scheme is less accurate than either the Brailovskaya or a Lax-Wendroff* differ- 
ence scheme. 

Victoria and Steiger (ref. 19) extended the original Crocco finite -difference scheme 
to two dimensions and solved the unsteady Navier-Stokes equations for the supersonic 
flow over a rearward -facing step. Good agreement was obtained between the numerical 
results and those found experimentally by Batt and Kubota (ref. 20). 


* A Lax-Wendroff difference technique is referred to here as an explicit method 
which has temporal and spatial truncation errors of second order. 
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MacCormack (ref. 21) has made calculations for the solution of a shock impinging 
on a laminar flat -plate boundary layer. In obtaining solutions to this problem, he modified 
a second-order finite -difference scheme which he had introduced earlier (see MacCormack 
(ref. 22)) for the problem of hypervelocity impact cratering. MacCormack's original 
method is a two-step Lax-Wendroff difference technique which alternately uses forward 
and backward differences for the convection terms. This scheme was modified in the 
incident -shock problem by splitting the governing equations into two sets - one for the 
x-derivatives and one for the y -derivatives. The advantage of the split system is that the 
computation proceeds with larger time increments since the stability criterion is less 
stringent. By using the split system, the required time step is the minimum time incre- 
ment of that found by applying the usual Courant -Friedrichs -Lewy conditions individually 
to the x and y subset of equations. 

MacCormack presented results for an incident oblique shock onto a flat boundary 
layer at M^ = 2.0 for which Hakkinen et al. (ref. 23) have made experimental measure- 
ments. The Reynolds number at the intersection of the shock wave and the flat plate, 
based on the length from the leading edge and free-stream conditions, was approximately 
3 x 10®. In general, the agreement between predicted and measured wall pressure and 
skin-friction distributions is good, although the x-grid spacing in the separation region 
is too coarse for adequate resolution. 
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C f,o 

c p,o 
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SYMBOLS 


constant of proportionality in linear viscosity law defined in equation (34) 
skin-friction coefficient, r 

skin -friction coefficient, 


1 

2 ' 


tP 0O U « 


ip_U^ 
2 M ° o 


pressure coefficient, 


P -Pp 
2 P ° u o 


speed of sound 

specific heat at constant pressure 
specific heat at constant volume 


E 


total energy 


G 

g 

k 

j,k 

L 

M 

N Pr 

P 

P 

P2 

Pt,2 

R 
R 

R 
R 

R 


internal energy 

amplification matrix 

eigenvalue of amplification matrix G 

thermal conductivity 

indices for x- and y -direction, respectively 

distance along flat plate from leading edge used in reference Reynolds number 

Mach number 

Prandtl number 

Curie correlation pressure 

pressure 

Stewartson and Williams’ correlation pressure 
total pressure behind a normal shock 
gas constant in equation of state 


R, 



Reynolds number, 

^ oo oo 

'A 

M oo 


Reynolds number, 

Poo U oo L 


Moo 


Reynolds number, 

P oo U =o X 

: °°,x 

Moo 

: °°> x c 

Reynolds number, 

Poo U =o X ( 

Moo 

•o,x 0 

Reynolds number, 

Po u o x o 

Mo 
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empirical constant in Sutherland viscosity law, for air S = 198.6° R 
temperature 

constant adiabatic wall temperature given in equation (38) 

free-stream stagnation temperature 

time 

time increment 

velocity component parallel to flat plate 
velocity components parallel and normal to ramp 
rarefaction parameter 
velocity component normal to flat plate 

skewed coordinate parallel to ramp; Curie correlation length 

Stewartson and Williams' correlation length 

Cartesian coordinates parallel and normal to flat plate 

distance between two successive grid points in x -direction 

skewed coordinate normal to flat plate 

distance between two successive grid points in y -direction 

ramp angle 

angle between dividing streamline and wall 
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y 

A 

6 

6 * 

M 

M B 

P 


<P 

Xoo 

Ip 


angle between u = 0 locus and wall 
ratio of specific heats, Cpy'Cy 
denotes either Ax or Ay 
boundary -layer thickness 

displacement thickness defined in equation (37) 

viscosity coefficient 

bulk viscosity coefficient 

density 

shear stress 

dependent variable 


weak -Interaction parameter, M 


nondimensional stream function 


3 \ I C °° 


JR 




Subscripts: 
c corner 

e edge of boundary layer 

o start of compression -corner interaction 

r reference 

w wall 


00 


free stream 



Superscripts: 


* temporarily used to denote nondimensional quantities 

n number of time cycles 


METHOD 

Governing Equations 

The governing equations that describe the motion of a viscous heat-conducting fluid 
are the Navier -Stokes equations which are given as follows for two-dimensional unsteady 
flow with respect to the Cartesian coordinates, x and y. 

Continuity: 


9£ a^u 9ev = o 

at ax ay 


x-momentum: 

8pu a(p + pu^) apuv _ 9t xx 9r xy 

at ax ay ax ay 


y -momentum: 


apy apuv 8 (p + pv 2 ) _ 9T yx 8r yy 

at ax ay ax ay 

Energy: 

aE au(E + p) av(E + p) _a_/ k ar\ _a_/, 8T\ 
at ax ay ax\ ax/ ay\ ay J 


* + VT xy) + + VT yy) 


where 


3u 


dX 




9u , 


b 3 A ax ay 


( 9u dv 

T Xy~Ty X - + — 


( 1 ) 


( 2 ) 


(3) 


(4) 


(Equations continued on next page) 
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In addition to these equations which express the conservation of mass, momentum, 
and energy, it is necessary to include a state equation which, for a perfect gas, is given 
by 

p = pRT (6) 

The viscosity coefficient is a function of temperature and is adequately approximated by 
Sutherland’s semiempirical equation 


p _ T r + S 0 / T \3/2 
Pr T +S 0 ^T r J 


where the constant S 0 has been experimentally determined to be 198.6° R for air. 

From equations (1) to (7) six equations in the six unknowns, u, v, T, p, p, 
and p are provided if the bulk viscosity coefficient Pg is set equal to zero. The 
effects of the bulk viscosity coefficient are important in sound propagation and shock- 
wave structure, as discussed by Vincenti and Kruger (ref. 24). In the present investiga- 
tion it is neglected since the grid spacing used in the region of the shock is too coarse to 
allow resolution of the shock structure. 


The variables in equations (1) to (7) may be nondimensionalized as follows: 

v 


x*=2L 

X L 


,*=L 


* - JL 

U 

oo 


V* = 


U 


3* = 

Poo 


p * = ^4 

P U " 

H oo oo 


T* = ^ 

" u 2 

OO 


E* = — — 


p u ■ 

“oo c 


II* = 




( 8 ) 


The resulting system of nondimensionalized differential equations may be conveniently 
written in vector form as 




+ 



(9) 
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where 


r p *^ 

. , P*u* 

w* = J 

)p*v* 




F* = 


r p*u* 
p* + p*u* 2 
\ p*u*v* 
u*(E* + p*) 


o* = 


"" p*v* 
p*u*v* 
p* + p*v* 2 

y^+p*). 


( 10 ) 


J 


and 


S* = 1 


R oo,L 


Np r 


9t, 


drZ 


xx “‘xy 


9x 


9y 


8T xy 8T yy 

ax* ay* 


d („* 

. * |p . * 

ax V ax 


a / ,* ar* 


®y 


m 


9y 


a / * * * * \ 

+ ^*( u T xx + V T xy) 


a / * * * * \ 

:^( U V +v T yy) 




T * =2 p* 8 ^-2 ^/V + *£\ 

XX ^ * O l 4c sk I 


ax* ^ \ax* ay* 1 


T* = n*/' 8u * | 

x y \ay* ax*; 


T * = 2u*^l 2p*/au* av*\ 

yy “ ay* " 3 U* ay*; J 


(11) 


( 12 ) 


Poo U oo L 
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with the Reynolds number L = — 77 and Prandtl number Np r = . In the pres- 

5 r" 00 K 

ent analysis, except as otherwise noted, the Prandtl number is maintained at a constant 
value of 0.72, which is the experimentally determined value for diatomic gases at moder- 
ate temperatures. The nondimensionalized Sutherland relation becomes 


P* = 


(y - i)m£ 

T* + S* 


+ £ T 


(y - Dm 2 t* 


2^l 3 / 2 


(13) 
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where is the free -stream Mach number and the constant S* is given by the 

equation 

s* = — (14) 

(y - 1 )M“T m 

It is observed that when the Sutherland viscosity is used, the dimensional free -stream 
static temperature is not eliminated from the governing equations and therefore must be 
specified in each calculation. The nondimensionalized state equation becomes 

p* - y ~ 1 p*T* (15) 

At this point the asterisk notation will be dropped with the understanding that all quantities 
are nondimens ional unless otherwise indicated. 


Finite -Difference Technique 

The finite -difference scheme chosen in the present investigation is the two-step 
explicit scheme proposed by Brailovskaya (ref. 2). Comparisons of steady -state solu- 
tions with solutions to Burger's equation (model of the Navier -Stokes equations) are pre- 
sented by Carter (ref. 18) between the Brailovskaya scheme, which has a truncation error 
of 0(At + Ax 2), and a Lax -Wendroff scheme with truncation error of 0(At2 + Ax2). Both 
schemes yield results of comparable accuracy near the asymptotic steady-state solution. 
Additional comparisons are presented by Carter (ref. 18) between solutions of the Navier - 
Stokes equations for a flat -plate flow field using the Brailovskaya scheme, and solutions 
obtained by Thommen (ref. 25) using a Lax -Wendroff scheme. Only small differences 
were found and hence it was concluded that these two schemes result in comparable accu- 
racy for obtaining steady solutions to the Navier-Stokes equations. The Brailovskaya 
scheme was chosen since the same grid is used for both time steps and therefore this 
scheme should be more efficient. In addition, variable grid was used for some of the 
calculations in the present investigation, and the Brailovskaya scheme is easier to modify 
to take noncentral differences into account. 

Application of the Brailovskaya scheme to equation (9) with t = n At, x = j Ax, 
and y = k Ay results in the following difference equations where the grid spacing has 
been assumed to be constant in the x- and y -directions: 


W ?S l " W j 1 ,k _ 

( F j+l,k " F j-l,k 

,-n r n 

G j,k+1 _G j,k-l 

At 

\ 2 Ax 

2 Ay 

w j l ,k 1 - W ij,k / 

V n+1 _F n+1 

F j+l,k F j-l,k 

n+1 — n+1 

G j,k+1 " G j,k-1 

At V 

] 

v 2 Ax 

2 Ay 


> (16) 
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In the first step, a temporary value of w ^denoted by is calculated at the new 

time step; this value is improved in the second step by reevaluating the convection term 
with the temporary values of w, the stress term S being repeated from the first step. 
The usual linearized stability analysis suggests that this two-step scheme is conditionally 
stable regardless of the magnitude of the Reynolds number. At the completion of each 
step the desired unknowns are computed from the vector w by the following equations: 



where the subscripts denote the elements of the column vector w. In equations (16) the 
viscous stress and heat conduction derivatives are contained in the term sf , . In the 
present investigation these terms were treated in their expanded form and approximated 
by central -difference quotients of second-order accuracy. For example, the shear term 
was approximated for constant grid spacing by 



An alternate procedure for this term is that used by Brailovskaya (ref. 2) and is written as 



(19) 

These approximations both result in the same order of truncation error; however, the 
expanded form (eq. (18)) is more efficient to use with variable grid and the skewed coor- 
dinate transformation which is used in the present investigation. 
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The Brailovskaya finite -difference scheme is conditionally stable since the maxi- 
mum time increment by which the solution may be advanced at any time step is dependent 
on the spatial grid size as well as the solution itself. An approximate analysis which 
yields an estimate of the maximum time increment is that of von Neumann which consists 
of examining the linearized difference equations for the amplification of short wavelength 
disturbances of small amplitude. Allen (ref. 26) presented a von Neumann stability 
analysis, based on a technique given by Richtmyer (ref. 27), for the Brailovskaya scheme 
applied to the inviscid equations. The resulting stability criterion was found to be the 
CFL (Courant -Friedrichs -Bewy) limit which is 
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This restriction on the value of At is a necessary condition for stability of the finite - 
difference solution of the Navier-Stokes equations if the calculations are made in inviscid 
regions.. Near the body surface where the viscous effects are predominant, a stability 
limit on At can be found which results from considering the viscous terms and is 
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The details of the analysis which yields these stability criteria are presented in appen- 
dix A. In performing the calculations, the minimum value of At given by equations (20) 
and (21) was used. 


Variable Grid 

Many flow fields contain regions which differ significantly in their characteristic 

lengths, and therefore difficulty is encountered in attempting to solve such flow fields 

with a uniform grid mesh. In the present investigation this particular problem arose in 

solving for the supersonic flow over a flat plate. As the distance x increases from the 

1/2 

leading edge, the boundary layer thickens like x ' , whereas the region bounded by the 
shock wave grows approximately like x. The net effect is that with increasing x, the 
viscous effects are confined to a smaller fraction of the distance from the wall to the shock 
wave. Since the variations in the inviscid shock layer are much less than those in the 
boundary layer, it is necessary to use a variable grid in order to attain maximum effi- 
ciency with the finite -difference mesh. 

A variable grid can be established either by transforming the independent variables 
with a suitable stretching function or by varying the grid in the physical plane. Skoglund 
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et al. (ref. 28) used a logarithmic coordinate transformation to open the grid in the region 
away from the intersection of an incident shock wave with a flat -plate boundary layer. 
Cebeci, Smith, and Mosinskis (ref. 29) expanded the grid at a constant rate away from the 
wall in their analysis of the turbulent boundary layer flowing over the body surface. In 
the present case, the latter technique of imposing a variable grid was used and the follow- 
ing difference quotients are easily derived from Taylor series along with their respective 
truncation errors: 
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Skewed Coordinate System 


Burstein (ref. 30) computed the inviscid supersonic flow in a channel which con- 
tained a compression ramp and utilized a Cartesian coordinate system throughout the flow 
field. The procedure resulted in extensive interpolation on the ramp since there the grid 



Figure 2.- Skewed coordinate system. 


points do not coincide with the ramp surface. This problem can be avoided by using a 
skewed coordinate system on the ramp as shown in figure 2. The skewed and Cartesian 
coordinate systems are related by 


X = x sec o> 

Y = y - x tan a 


(26) 


The equations relating the derivatives with respect to the Cartesian coordinates (x,y) and 
those of the skewed coordinates (X,Y) are 
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As shown in figure 2, the Cartesian and skewed coordinates are used along the flat plate 
and ramp, respectively; however, special consideration has to be given to the interface 
between the two systems located at x = X = 0. In particular, the derivatives with respect 
to x (or X) require special treatment in order to maintain second-order accuracy when 
the difference equations are formulated along the interface. These difference expressions 
are obtained in the usual manner from Taylor series expansions and are presented in 
appendix B. 


Flat -Plate Boundary Conditions 

Figure 3 shows a schematic diagram of the flat -plate flow field at the leading edge. 
The flow field to be computed is enclosed in a box as shown in figure 3. The free -stream 
conditions are specified along the upstream boundary and also the outer boundary, pro- 
vided it is located outside of the leading-edge shock wave. Along the downstream bound- 
ary, the flow variables are unknown and must be evaluated from the oncoming flow. In 
the present calculations, quadratic extrapolation of the flow variables in the x-direction 
was used continuously to update the downstream boundary conditions. The wall boundary 
conditions are specified as shown in figure 3. It is observed that the wall pressure is 
unknown and must be evaluated from the neighboring flow field. Since the wall forms a 


Free -stream 
conditions 



Figure 3.- Schematic diagram of the flat plate flow field 
and computational boundaries. 
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boundary of the grid system, the central difference formulation of the equations cannot 
be used there to find the pressure as is done away from the wall. Different methods of 
computing the wall pressure and numerical tests of the downstream boundary conditions 
are discussed in the section "Results and Discussion." 


An obvious simplification in these boundary conditions is neglect of the velocity slip 
and temperature jump that occur at the wall near the leading edge. These effects are 
important in the merged-layer region (shock wave not distinct from boundary layer) which 
extends from the leading edge downstream to ~ 0.15 where 


V 


CO 


= M„ 



(28) 


according to measurements made by McCroskey, Bogdonoff, and McDougall (ref. 31). 
Downstream of the merged layer for supersonic Mach numbers is the weak-interaction 
region in which the effects of slip and temperature jump are negligible. Because of the 
higher Reynolds number range, it is the weak-interaction region that is of interest in this 
investigation; hence, for simplicity, the details of the rarefaction in the relatively small 
merged-layer region are assumed to have only a local effect. This assumption was 
verified a posteriori by comparing the numerical solution with results obtained from 
weak-interaction theory and similar solutions to the boundary -layer equations. 


Compression-Corner Boundary Conditions 

Figure 4 shows a schematic diagram of the computational box used in the 
compression -corner calculations with the boundary conditions indicated on the respec- 
tive faces of the box. These boundary conditions are the same as those used in the flat- 
plate calculations except for those on the outer boundary. In order to reduce the com- 
puter requirements of the calculations, the outer boundary was placed between the wall 
and the leading -edge shock wave. Placing the outer boundary in the disturbed part of the 
flow field requires that extrapolation be used continuously to update the flow variables 
since the conditions are unknown along this boundary. For the present calculations, the 
extrapolation on the outer boundary was based on the approximate simple -wave character 
of the outer inviscid flow field. The procedure used in this extrapolation can be described 
in the following manner. At the completion of each step of the two-step Brailovskaya dif- 
ference scheme, the inclinations of the steady left-running characteristics were computed 
along the first row of grid points inside the outer boundary. If the flow is assumed to be 
of the simple -wave type, these characteristics were linearly extended to the top row of 
grid points. The flow properties are invariant along straight characteristics; therefore, 
the unknowns at the top grid points were found from linear interpolation of those values 
assumed to be constant on the characteristics at the next to top row. 




Figure Schematic diagram showing the computational box and boundary 
conditions for the compression-corner calculations. 

The flow field along the outer boundary is of the simple -wave type provided that 

(a) The strength of the reflected waves from both the leading -edge shock and the 
vorticity layers generated by the leading-edge shock is negligible, and 

(b) The coalescence of the compression waves from the corner flow occurs beyond 
the outer boundary. 

It is reasonable to assume that the first condition is satisfied based on the success of the 
shock-expansion technique, which ignores these effects. Waldman and Probstein (ref. 32) 
showed that the strength of a reflected wave at a shock is less than 1 percent of the inci- 
dent wave for y - 1.4, free -stream Mach numbers less than 4, and flow deflection angles 
less than 10°. Even for an infinite Mach number, the strength increases to only 14 per- 
cent for deflection angles less than 44°. The vorticity generated by the leading -edge 
shock is negligible since its curvature is very small for the present conditions. Satisfac- 
tion of the second condition had to be verified a posteriori, although one would expect the 
shock formation point to occur at distances greater than two or three boundary -layer 
thicknesses, which was the typical position of the outer boundary. Clearly, this assump- 
tion is limited by Mach number since as the free-stream Mach number increases, the 
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separation and reattachment shocks lie closer to the surface, as has been discussed by 
Holden (ref. 33). 


RESULTS AND DISCUSSION 


Moo = 3.0 Flat -Plate Calculations 


Calculations were made for the supersonic, viscous flow over a flat plate with the 
following free -stream conditions: 

Moo = 3 -° R oo,L = 1q3 N Pr = °* 72 y = 1A 

The Sutherland viscosity law was used with the dimensional free -stream temperature 
chosen at 390° R. The initial conditions were chosen to be that of a flat plate impulsively 
accelerated to free -stream conditions. Free -stream boundary conditions were enforced 
along x = 0 (with the exception of the grid point located at y = 0 where the wall condi- 
tions were imposed) and along the outer boundary as shown in figure 3. The boundary 
conditions along the wall y = 0 are given by 


u(x,0) = v(x,0) = 0 

(y - 1 )m; 


T <*. 0) = T o,«, = 7— 1^2 (‘ * V M ”) 




(29) 


an isothermal walk being assumed. The density at the wall was found by quadratic 
extrapolation in the y-direction according to the relation 


p(x,y) = Pw( x ) + Pj(x) y + p 2 (x) y 2 


with the quantities p w , Pp and pg evaluated with the latest known values of the density 
at the points y = Ay, 2 Ay, and 3 Ay at each x-station. Other methods of finding the 
wall density or pressure (only one of these two quantities is needed for the isothermal 
condition since the other is determined from the state equation) are discussed later in 
this section. The values of the dependent variables along the downstream boundary, 
which was placed at x/L = 1.5, were continuously updated by quadratic extrapolation in 
the x-direction. Numerical tests on the effect of this approximation are also discussed 
later in this section. 


Computed results for different grid sizes .- The calculations outlined were made 
for three different grid sizes in order to obtain a measure of the required resolution. 
For the isothermal condition it was found necessary to have equal grid spacing in the x- 
and y -directions in the immediate vicinity of the leading edge so that numerical instabil- 
ity would not result in a divergent solution. It should be noted that there are inconsis- 
tencies in the literature with respect to this point. For similar calculations, Kurzrock 
(ref. 34) observed the same result as found here, whereas MacCormack (ref. 21), in 
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computing the supersonic flow over an adiabatic sharp-edged plate, used a grid spacing 
ratio Ax = 40 Ay and did not encounter numerical instability. Thommen's calculations 
(ref. 25) were also for an adiabatic flat plate, but unlike MacCormack, he observed that 
for Ax = 2 Ay, an instability resulted that was eliminated for Ax = Ay. 


In the present calculations the grid spacing was set equal to 0.05, 0.025, and 0.015. 
In the first two cases the grid was maintained constant throughout the computational box, 
but in the third case only Ay/L was held constant at 0.015, whereas Ax/L was varied 
as follows in order to reduce the number of x grid points from 100 to 66: 


0.015 

Li 


Ax 


= 0.020 


^ = 0.025 
Li 


(osf.so.15)' 
(o.l 5 s 0.25) 
(0.25 ^ s 1.5) 


( 30 ) 


The converged results for the three grid sizes are presented in figures 5,6, and 7. 
The effects of the coarseness of the grid and incorrect wall boundary conditions in the 
leading -edge region are clearly shown in figure 5 by the large oscillations in the wall 
pressure near the leading edge. These oscillations disappear closest to the leading edge 
for the 0.015-grid solution; further downstream, the oscillations disappear in the 0.025- 
grid solution as these results approach the 0.015-grid results. The results obtained with 
0.05 grid differ significantly from the finer grid results. For the leading -edge calcula- 
tions the parameter which determines the effective resolution of the calculations is the 
grid-spacing Reynolds number = P °°^°° — . These calculations were made for 

R^ = 15, 25, and 50 for the three respective grid sizes. Numerical instability resulted 
when the same flow field was calculated with R^ = 250. Therefore, for the leading-edge 
calculation, there appears to be an upper bound on the grid-spacing Reynolds number in 
order to achieve stable results. 

Figure 7 shows profiles of u/u^, v/u B , T/T m , and p/p rj0 at x/L = 1.0 
( R oo,x = 10 ^). The extrapolation at the downstream boundary x/L = 1.5 introduces 
slight errors, and therefore it was desired to make the comparisons upstream of the 
influence of this effect. In figure 7 the distinctive features of the flow are evident. The 
boundary -layer edge is at y/L ~ 0.25 with the leading -edge shock located at y/L ~ 0.57; 
thus, x/L =1.0 is downstream of the merged region. This result is expected since 
V M = 0.083 at this position and, as indicated previously, the downstream extent of the 
merged region is ~ 0.15. In figure 7(c) the slightly negative wall temperature gra- 
dient indicates that heat is being transferred from the wall to the stream. This result is 
expected since the adiabatic wall temperature for Np r = 0.72 is less than the free- 
stream total temperature which was the isothermal wall condition chosen here. 
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Convergence of the calculated results with decreasing grid size is evident in the 
profiles at x/L = 1.0 shown in figure 7. Observing the boundary -layer part of these 
profiles, it appears that it is necessary to use approximately 15 grid points in the bound- 
ary layer for adequate resolution. This result should be general and is used as a guide- 
line throughout the present investigation. Cheng (ref. 14) determined qualitatively from 
truncation-error considerations of Burger’s equation that for second-order accurate 
schemes, a minimum of 20 or 30 grid points are required per characteristic viscous 
dimension for adequate resolution. Thommen and Magnus (ref. 35) concluded from their 
numerical study of Burger's equation that only 5 grid points would be required per 
boundary -layer thickness. 

The effect of the grid refinement is more significant near the shock wave as shown 
in figures 7(b) and 7(d) for the profiles of normal velocity and density. The 0.015-grid 
spacing eliminates most of the oscillation in the shock transition region, whereas the 
0.025-grid size results in significant oscillation. Naturally, this shock smearing places 
some upper bound on the grid spacing; however, as the calculations are extended down- 
stream, this grid criterion can be relaxed somewhat since the leading -edge shock weakens, 
and the shock layer thickens and allows more distance for the oscillations to damp out 
before entering the boundary layer. 

Numerical results 
Grid size 
O .050 
□ .025 
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Figure 5* - Comparison of wall pressure distribution for different 
grid sizes with weak interaction theory. 
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Comparison with weak-interaction th eory.- Shown for comparison in figure 5 is the 
wall pressure predicted by the second-order weak-interaction analyses of Lees and 
Probstein (ref. 36) and Kubota and Ko (ref. 37). The analysis of Lees and Probstein is 
based on a Taylor series expansion of the induced pressure in powers of d6*/dx with 
the coefficients in this expansion determined from the tangent -wedge approximation. 

The weak-interaction analysis of Kubota and Ko is somewhat similar in that it is based 
on a series expansion in powers of where 

with the coefficients in the series determined by substitution into the integral boundary- 
layer equations of Lees and Reeves (ref. 4). In addition, the equations given by Kubota 
and Ko (ref. 37) were for Np r = 1.0 and an adiabatic wall. With Np r =1.0 the 
approximate recovery factor is 1.0 and results in an adiabatic wall temperature equal to 
the free -stream total temperature, the assumed wall temperature in the present calcula- 
tions. As seen in figure 5, both theories give good results in comparison with the present 
numerical results. The Lees and Probstein (ref. 36) theory slightly underpredicts the 
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Figure 6.- Comparison of wall skin-friction distribution for 
different grid sizes with weak interaction theory. 


.4 



induced pressure which is consistent with the comparisons made by Kubota and Ko 
(ref. 37). However, a much greater disagreement is shown in figure 6 where the wall 
skin -friction coefficient written as 


c f,w = 


Tw 



( 32 ) 


is plotted against x/L and the interaction parameter Xoo* Excellent agreement is 
obtained between the Lees and Probstein analysis and the numerical results, but the 
agreement is much poorer between these results and the Kubota and Ko analysis. The 
difference is unexplained, although downstream the agreement improves considerably, as 
shown by Carter (ref. 18). 


C omparison with similar solutions .- Shown for comparison in figure 7 are the zero- 
pressure -gradient similar solutions of the boundary -layer equations for the specific flow 
conditions prescribed here. In order to convert the similar -solution results given in 
tabular form by Low (ref. 38) from the nondimensional Blasius variables to the present 
nondimensional variables, it is necessary to specify u, p, and T at the boundary -layer 
edge as well as u = v = 0 and T = T w at the wall. Rather than choose free -stream 
conditions for the outer-edge conditions as would be consistent with first-order boundary- 
layer theory, the edge conditions were chosen from the present results obtained from the 
Navier-Stokes equations. This choice was made by assuming that the wall pressure ratio 
Pw/P°° = was constant throughout the boundary layer. This pressure and the edge 

temperature Tg/T^ = 1.171, taken from figure 7(c), gives the edge density Pq/P^ = 1.221. 
From figure 7(a) the edge velocity is Ug/u^ = 0.95. This procedure obviously forces the 
agreement of the Navier-Stokes and boundary-layer equation solutions at the boundary - 
layer edge; however, the good agreement elsewhere indicates the correctness of the pres- 
ent results. The effect of the negative pressure gradient is seen in figure 7(a) since the 
u-component of velocity is greater at a given y-position for the Navier-Stokes results than 
that for the zero-pressure-gradient similar solution. In addition to the differences due to 
pressure gradient, it was also found that the Chapman -Rubesin (ref. 39) form of the linear 
viscosity law used in the similar solutions introduced a slight error. This viscosity law 
is given by 



(33) 


where 


C = 


^r + So 

T w + S 0 



(34) 


At the wall the linear result is forced to agree with that given by the Sutherland relation. 
Inserting the Chapman-Rubesin viscosity law into the Navier-Stokes calculations reduced 
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the differences between the Navier -Stokes results and the similar solutions by several 
percent. 

Computation rate . - It is convenient at this point to discuss the computer time used 
to perform the present calculations. The large demands for both computer time and 
storage which are required for these calculations are the main drawback to the present 
approach. The computational rate of the present program as applied to the flat -plate 
flow field was found to be 2.75 X 10® grid points/hour on the Control Data 6600 computer 
at the Langley Research Center. At this rate the program can execute 10® time cycles 
through a field of 2750 grid points in 1 hour of machine time. 

The number of cycles required for convergence of the calculations discussed were 
approximately 500, 900, and 1300 for the grid sizes of 0.05, 0.025, and 0.015. In all three 
cases the initial conditions were assumed to be those of a flat plate impulsively accel- 
erated to free-stream conditions. The required computer times for these calculations 
for the respective grid sizes were 0.175, 0.82, and 2.12 hours on the Control Data 6600 
computer. For the flat -plate calculations, the two velocity components at the downstream 
end of the calculation were the slowest of the dependent variables to converge. Conver- 
gence was assumed when their values ceased to change in the fifth significant digit. As 
is typical of results obtained by the time -dependent technique, the solution varies rapidly 
at the start, and is followed by a slowly varying monotonic approach to convergence. 

Different methods of computing wall pressure . - Several techniques for updating the 
fluid density (or pressure) at the wall were tested on the leading-edge flat -plate flow field 
discussed previously. This calculation provides a severe test of the various numerical 
methods because of the large gradients near the leading edge. 

The techniques that were used fall into two categories: extrapolation from interior 
points, and evaluation of the equations at the wall with one-sided difference quotients. As 
was indicated previously, quadratic extrapolation of the density from the three grid points 
above the wall and the subsequent evaluation of the pressure from the state equation gave 
the best results. Linear extrapolation of the density resulted in an unstable calculation. 
This result is not surprising since linear extrapolation introduces an error of 0 (a 2) in 
the dependent variable, whereas the difference scheme computes the dependent variable 
with a higher order error of 0( A®), as can be seen in equation (16). Quadratic extrapola- 
tion admits a third-order error and hence is consistent with the difference scheme. In 
addition to the density extrapolation, quadratic extrapolation was also attempted with p 
and with p + pv^ to find the wall pressure. Both of these attempts produced unstable 
calculations. 

In principle, it would seem to be preferable to use either the continuity or 
y-momentum equation to determine either the density or pressure at the wall. Evaluated 
at the wall, the continuity equation becomes 
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( 35 ) 


dp _ 8py 

at ay 

Similarly , the y -momentum equation is given by 
fe + pv 2 _ 1 

9y R oo,L 


Ji A + IJjl. 
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(36) 


Equation (35) or (36) is evaluated to find p or p at the wall, respectively, by using 
suitable second -order one-sided difference quotients for the spatial derivatives and a 
first-order forward time -difference quotient for 3p/9t. The calculations were unstable 
when the y -momentum equation was used to evaluate p w ; for the same calculations, use 
of the continuity equation to find p w resulted in a converged solution although large 
oscillations occurred near the wall. Figure 8 shows a comparison of the pressure dis- 
tributions at x/L = 1.0 obtained by using extrapolation and the continuity equation to 
find p w . In both calculations the grid spacing was ~ = = 0.025. Using the con- 

tinuity equation to find p,„ results in a pressure distribution with large oscillations 
whose amplitude decreases away from the wall. Clearly, the extrapolation technique is 
preferable. 

Numerical tests on downstream boundary conditions .- The use of extrapolation in x 
continuously to update the flow variables at the downstream boundary presumably approx- 
imates the effect of a flat plate which extends far downstream. Naturally, the questions 
arise: What is the degree of inaccuracy introduced by this approximation and what is the 
extent of its upstream influence through the subsonic part of the boundary layer? It is 
desired to extend the present calculations downstream to higher Reynolds numbers by 
using the present downstream conditions as upstream conditions for the next calculation. 

If the extrapolation introduces sizable upstream errors, then such a procedure would not 
be efficient because of the large degree of overlapping required for the computational 
boxes. Callens (ref. 40) discussed the use of such a procedure and suggested that over- 
lapping would reduce the effect of the errors due to extrapolation; however, he did not 
indicate the magnitude or upstream extent of these errors. 

A numerical test was performed on the calculations discussed previously in this 
section to answer these questions. Calculations were made by use of the 0.025 grid for 
three computational boxes, one of which enclosed the other two. The results of these 
calculations are shown in figure 9, where the x-extent of each calculation is given. Both 
the wall pressure and the pressure along y/L = 0.225, which is above the sonic line and 
near the boundary -layer edge, are plotted against x/L in order to compare the effect 
of the error on both a subsonic and supersonic region. The "bump" in the pressure dis- 
tribution along y/L = 0.225 is part of the oscillations that result from the smearing of 
the shock wave by the finite -difference scheme. The results shown in figure 9 indicate 
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Figure 9-- Effects of downstream extrapolation on streamwise 
pressure distribution near the downstream boundary. 

that an error of about 1 percent in the wall pressure occurs at x/L = 1.0 where the 
quadratic extrapolation was made. This error is propagated upstream through the sub- 
sonic part of the boundary layer and disappears within approximately 10 grid points. The 
effect of the extrapolation on the pressure along y/L = 0.225 is considerably smaller 
although a slight error is detectable. This error is attributed to the fact that the point 
x/L = 1.0, y/L = 0.225 is slightly within the zone of influence of the boundary layer 
which is altered by extrapolation errors. 

The second aspect of the test calculations is to determine the effect of using 
extrapolated conditions as upstream boundary conditions for a downstream flat -plate cal- 
culation. Figure 9 shows that the wall pressure quickly recovers to the correct solution, 
whereas the slight error in the upstream pressure along y/L = 0.225 persists even 
farther downstream as is typical of inviscid flow. 

The important point to be made here is that the use of extrapolation introduces only 
small errors of limited upstream extent. It is presumed that as the calculations extend 
farther downstream, these errors will become less since the gradients in the x -direction 
become smaller. In addition, it does not appear to be necessary to overlap the computa- 
tional boxes since the extrapolation errors are quickly damped out. 

Extension to higher Reynolds numbers .- The present results were extended down- 
stream to higher Reynolds numbers by using the downstream conditions at x/L = 1.0 as 
upstream boundary conditions as was discussed in the previous section. In this manner 


33 


I 


the calculations can be made for any Reynolds number desired. Instead of using a rec- 
tangular computational box as was done in the leading -edge calculations, a trapezoidal 
box was used with the outer boundary of the downstream box inclined at an angle slightly 
larger than the leading-edge shock-wave angle at the interface. The expansion waves 
generated by the displacement thickness weaken the shock, and its position therefore 
moves farther inside the computational box. Use of the trapezoidal box avoids the need- 
less calculations in the free stream that are required for the rectangular box and also 
allows free -stream conditions to be specified along the outer boundary. 

The calculations were made from x/L = 1.0 (k^x = 10®) to x/L = 6.0 
x = 6.0 x 10®) by use of the trapezoidal computational box with the outer boundary 
inclined at an angle of 29° with respect to the free stream. The grid spacing in the 
y-direction along the inclined boundary was increased at a constant rate of 20 percent 
from 0.025 to 0.055. Since the grid spacing in the x-direction is determined by the 
y-spacing along the inclined boundary, the resulting x-grid increased from 0.025 to 0.1. 
Figure 10 shows a flow -field map deduced from the leading -edge and downstream flat- 
plate calculations. The characteristic features of the supersonic flat -plate flow field 
are shown: the boundary -layer -induced shock wave, displacement thickness, streamline 
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Figure 10.- Computed flow field over a flat plate with 
M m = 3*0 and T w = To,°°* 
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pattern, and contours of constant pressure. The displacement thickness was computed 
by using the trapezoidal rule to evaluate 

6 *=i'ot 1 -<^] dy <37) 

It was observed that the u -velocity component reaches a local maximum near the boundary - 
layer edge, as can be seen in figure 7(a). By denoting this value as u e , 6 was taken as 
the value of y where u = 0.99u e . Also shown in figure 10 is the excellent agreement 
obtained between the calculated displacement thickness and that given by the Kubota and Ko 
weak-interaction theory applied to the present flow field. The finite -difference scheme 
smears the shock over several grid spaces (the actual number depends on the grid size), 
and therefore the discrete shock wave shown in figure 10 is an estimate deduced from the 
pressure distributions. It was found that the scatter in the locus of estimated shock 
points was always no more than one x- or y-grid spacing. The pressure contours cor- 
rectly indicate both the relatively constant nature of the pressure in the boundary layer 
and the constant slope that is characteristic of simple -wave flow in the inviscid part of 
the flow field. The inclination of the pressure contours differ from that of the Mach lines 
by less than 0.5° in the inviscid part of the flow exclusive of the immediate neighborhood 
of the shock. 


= 6.06 Flat -Plate Calculations 

OO 

In the compression-corner calculations, which will be described in the next section, 
comparison is made with the experimental measurements of Lewis, Kubota, and Lees 
(ref. 11). The upstream boundary conditions for this calculation were computed by use 
of five tandem computational boxes which extend from the leading edge to x/L = 10.0 
(«»,* = 105), as shown in figure 11. Greater use was made of variable grid in these cal- 
culations than in the Moo = 3.0 calculations, as shown in the table given in figure 11. 
Variable grid in both the x- and y-directions was imposed by varying the respective grid 
spacing at a constant rate, the rate being 10 percent or less. The effect of variable grid 
is easily seen in figure 11 since the first four boxes have approximately the same number 
of grid points, but their respective sizes increase considerably with downstream distance. 
As seen in figure 11, no overlap of the computational boxes was used as was previously 
found to be acceptable. 

For these calculations, M m = 6.06 and the free -stream static temperature was 
88° R. The wall -temperature boundary condition was assumed to be best approximated 
by an isothermal condition since longitudinal conduction occurred in the experimental 
model. (See ref. 41.) If the recovery factor is assumed to be equal to \/Np r , it is found 
from the oblique shock equations that the recovery temperature only changes by 2 percent 
from the flat plate to the ramp; thereby, the use of an isothermal condition is further 
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Figure 11.- Schematic diagram of the computational boxes used to 
compute the = 6.06 flat -plate flow field. 

justified. The constant wall temperature is denoted by T aw and is given by 
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The results obtained here are similar in nature to those computed for = 3.0 
and therefore only a short discussion will be given. Figure 12 gives a flow-field map 
showing the shock wave, streamlines, pressure contours, and displacement thickness. It 
is seen in figure 12 that the shock position found numerically is in excellent agreement 
with that found by Lewis et al. (ref. 41) from pitot pressure measurements. Also shown 
is the displacement thickness predicted by the Kubota and Ko weak-interaction analysis, 
which is slightly larger than the computed displacement thickness. Correspondingly, the 
estimate of the induced pressure based on the Kubota and Ko analysis is slightly larger 
than that found in the present investigation, as shown in figure 13. It is presumed that 
both of these differences are due to the fact that the Kubota and Ko analysis assumed that 
Np r = 1. This value yields a higher recovery temperature that results in a thicker 
boundary layer and greater induced pressure. In figures 13 and 14, it is observed that 
the induced pressure and the skin friction found by solving the Navier -Stokes equations 
has slight oscillations near x/L = 1.6 and x/L = 4.75. These oscillations are non- 
physical and result from the use of extrapolation at the downstream boundary as was 
previously discussed. Also shown in figure 13 is the wall pressure measured by Lewis 
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Figure 12-.- Computed flow field over a flat plate with 
= 6.06 and T w = Taw* 

et al. at R^x = 7.5 X 10^ and it is seen to be 7 percent higher than that found numeri- 
cally or that given by weak-interaction theory. 

Figure 14 shows the wall skin -friction coefficient plotted against x/L and the 
viscous interaction parameter Xoo- The agreement between the Lees and Probstein weak- 
interaction analysis and the present results is excellent, whereas the Kubota and Ko anal- 
ysis overpredicts the increase in skin friction, as was found in the M«> = 3.0 compari- 
sons. Also shown for comparison is the similar solution skin-friction distribution, which 
is less than the other results since ignoring the favorable pressure gradient results in a 
thicker boundary layer and reduced skin friction. 

Figure 15 shows a comparison of pitot pressure from the present results with that 
measured by Lewis et al. at R M x = 7.5 x 10^, and it is seen that the agreement is only 
fair. Inadequate numerical resolution is not the cause of the differences shown in fig- 
ure 15 since the calculated results changed by less than 1 percent when the grid spacing 
was halved in the direction normal to the plate. 

Figure 16 shows a similar comparison of the numerical and experimental Mach num- 
ber profile at the same x-station. The Mach number has been normalized by the edge 
value which was determined to be 5.56 in the experiment and 5.66 from the numerical 
solution. Similarly, the distance normal to the plate has been nondimensionalized by the 
boundary -layer thickness which is defined in figure 15 and is the same definition as that 
used by Lewis et al. The numerical and experimental values of the boundary -layer 


37 










Present results 



thickness are 2.31 mm (0.095 inch) and 2.18 mm (0.086 inch), respectively. The agree- 
ment of the two profiles shown in figure 16 is excellent. Plotting the profile data in this 
manner enhances the agreement, but it should be noted that the improved agreement in 
figure 16 is better than that shown in figure 15 because M/M e is less sensitive than 
p t 2 ^P 0 ^ For Mach numbers near 5, a 1 -percent error in Mach number amplifies to a 
2 -percent error in p t 2^0 «>' 

Me© = 3.0 Compression-Corner Calculations 

The Moo = 3.0 flat -plate calculations which were made downstream to 
Hoo,x = 6.0 x 10 3 were used as upstream boundary conditions for the compression-corner 
calculations with the corner located so that the Reynolds number based on free-stream con- 
ditions and the distance from the leading edge to the corner is Rc©,x c = 1-68 x 10^. The 
outer position of the computational box was located at y/x c = 0.0868 which is approxi- 
mately one -half of the shock-layer thickness and twice the boundary -layer thickness at 
the upstream boundary. The computational box extended downstream to x/x c = 1.99. 

The total number of grid points used was 2156 (77 in the x-direction and 28 in the 
y-direction); thereby a constant grid spacing in the x- and y-directions of 0.0214 and 
0.00321, respectively, resulted. For the compression-corner calculations, the computa- 
tional rate was 2.25 X 10® grid points/hour on the Control Data 6600 computer. This rate 
was a decrease of approximately 20 percent from that found for the flat -plate calculations 
and was due to the extra calculations needed at the coordinate system interface and for 
the simple -wave extrapolation. The same convergence criterion was applied that was 
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Figure 16. - Comparison of theoretical and experimental Mach number profile. 

used in the flat -plate calculations. The compression-corner calculations to be discussed 
in the remainder of this section typically required 1500 to 3000 time cycles to converge, 
depending on the initial conditions. 

The resulting wall pressure distributions for T w = T 0 and a = 5.0°, 7.5°, and 
10.0° are shown in figure 17, where the flat -plate distribution is based on the Lees and 
Probstein weak-interaction theory. It is seen that the pressure continues to decrease for 
some distance from the upstream position of the computational box; thereby, it is indicated 
that all the upstream influence of the compression corner is contained. Flow separation 
occurred in the 10.0° case; and the points of separation and reattachment are denoted by 
S and R, respectively. It is observed that two more inflection points occur in the wall 
pressure distribution in the separated case than in the unseparated cases. These inflec- 
tion points are the onset of the plateau region which occurs when the extent of separation 
is greater. The overshoot of the pressure above the inviscid level is expected since 
further downstream the static pressure will decrease because of the influx of streamlines 
having lower stagnation pressure due to shock losses. 
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Figure 17*- Comparison of wall pressure distributions for 
ramp angles of 5-0°, 7*5°j and 10.0°. 

Figure 18 shows the variation of the displacement thickness measured normal to 
the surface. The same definition of displacement thickness that was used in the flat- 
plate calculations is used here; however, the edge quantities are taken as those just inside 
of the outer boundary of the computational box. Plotting the displacement thickness in 
this manner gives a clear indication of the effect of the compression process on the bound- 
ary layer. Near the upstream boundary, the displacement thickness is seen to be in good 
agreement with that predicted by Kubota and Ko weak-interaction theory for a flat plate. 
Farther downstream, the upstream effects of the corner result in a rapid increase in the 
displacement thickness with it reaching a maximum at the corner. Past the corner, the 
boundary -layer edge approaches the surface until a minimum thickness is reached at the 
neck, after which the boundary layer continues to grow in a weak-interaction flow field at 
a new Mach number. It should be noted that at the corner, the displacement thickness was 
computed with respect to both surfaces, but the two values did not differ to within the 
plotting accuracy of figure 18. 
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Figure l8.- Comparison of the displacement thickness measured normal to 
the surface for ramp angles of 5 - 0 °> 7*5°j and- 10.0°. 

In figure 19, the wall skin -friction distributions for the ramp angles of 5.0°, 7.5°, 
and 10.0° are compared. At the corner, the skin friction is double -valued since at this 
point the shear stress is considered with respect to two surfaces. Also shown is the 
skin-friction coefficient predicted by the Lees and Probstein weak -interaction theory for a 
flat plate. The skin -friction distribution clearly shows the extent of the separation region 
for the 10.0° case. It is also seen that the 7.5° case for these flow conditions is that of 
incipient separation since the skin -friction coefficient is very close to zero at the cor- 
ner. Comparison of figures 18 and 19 shows that the flat -plate weak-interaction estimate 
of the downstream skin -friction coefficient agrees with the computed results, whereas 
the estimated displacement thickness does not. This result is expected since, to first 
order, the skin-friction coefficient depends only on the square root of the Reynolds num- 
ber, whereas the displacement thickness depends on both the square root of the Reynolds 
number and the square of the Mach number. 

A computed flow -fie Id map in the immediate region of the corner is shown in fig- 
ure 20 for the a = 10.0° case. The streamline pattern is shown inside as well as exte- 
rior to the separation bubble. There are approximately 17 grid points along the wall from 
separation to reattachment; thereby adequate resolution is given. In the separation bub- 
ble, the locus of u = 0 is shown to indicate the regions of forward and reverse flow. 
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Figure 19- - Comparison of wall skin-friction distributions for 
ramp angles of 5-0°, 7.5°, and 10.0°. 

Comparison of the values of the stream function indicate the low mass flow in the recir- 
culation region. It is also observed that the points of separation and reattachment are 
not equidistant from the corner as the bubble extends farther downstream from the corner 
than it does upstream. 

In the flow exterior to the recirculation region, the computed displacement thick- 
ness is plotted and is seen to be nearly parallel to the neighboring streamlines, which is 
consistent with the displacement body concept long used in boundary -layer theory. The 
pressure contours shown in figure 20 indicate that the pressure variations are small 
between the wall and the displacement thickness but larger between the displacement 
thickness and the boundary -layer edge. Outside the displacement thickness, the pressure 
contours quickly turn to follow the compression waves generated by the turning process. 
The pressure contours converge only slightly from the wall to the outer boundary of the 
computational box and hence the shock formation occurs exterior to the region of 
calculation. 

Effect of suction .- A calculation was made for the M m =3.0 and a = 10.0° case 
in order to determine the effect of wall suction on the compression-corner flow. Suction 
reduces the extent of the separated region by removing the section of the boundary layer 
with low total pressure; thus the boundary layer is thinned and the skin friction is 
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Figure 20 .- Computed flow field over a 10 ° compression corner for = 5-0- 

increased. The wall suction velocity was set at v/uoo = -0.01 for 0.786 = £- = 1.214, 

A c 

which is a region comparable to the length of separated flow 0.84 =£-= 1.22, which 

x c 

was found for the same case with no suction. 

The effects of the suction boundary condition are very evident in comparing the 
flow-field maps in figures 20 and 21. As a result of the mass removal, the pressure rise 
and streamline turning generated by the compression corner occur much closer to the 
corner. Examination of the skin-friction distribution shows that this arbitrary choice of 
suction velocity results in a case of incipient separation since Cf >w = 0 only at the cor- 
ner. Comparison of the displacement thickness between figures 20 and 21 shows the 
boundary -layer thinning produced by the suction. The kinks in the streamlines at the 
initial station, where the suction is first applied, are nonphysical and are a numerical 
result of the discontinuity in the boundary conditions. It was found that with the 1 -percent 
wall suction velocity, the amount of fluid removed in the corner region was 3.22 percent 
of that flowing through the upstream boundary. 

Numerical tests of simple -wave extrapolation .- Several numerical tests were made 
to determine the errors introduced into the calculations by the simple -wave extrapolation 
that was used to update the boundary conditions along the outer boundary. First the tests 
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Figure 21.- Computed, flow field over a 10° compression corner with 
a discontinuous wall- suction velocity of v w /Uoo = -0.01. 

were performed on a known flow field, namely, the Moo = 3.0 flow over a flat plate shown 
in figure 10. Two different positions of the outer boundary, y/L = 0.75 and y/L = 1.3, 
were used, and the resulting pressure distributions along these lines are compared in 
figure 22 with the solution obtained with the outer boundary located in the free stream. 

In both cases, the calculations were extended upstream to x/L = 2.0 and downstream to 
x/L = 5.0. Several observations can be made about the results shown in figure 22. First, 
it is seen that the calculation with the outer boundary at y/L = 0.75 results in a pres- 
sure distribution (which is the extrapolated pressure) that is in excellent agreement with 
that given by the calculation with the exact outer boundary conditions. Comparison of the 
pressure distribution along y/L = 1.3 given by the simple -wave extrapolation with that 
found by using the exact outer boundary conditions indicates differences up to 5 percent. 

It is also seen in figure 22 that the errors introduced along y/L = 1.3 by the simple - 
wave extrapolation reflect back into the flow field as erroneous expansion waves as 
deduced from the pressure distribution given by this computational box along y/L = 0.75. 
The first point of influence of the outer boundary at y/L =1.3 on the pressure along 
y/L = 0.75 is approximately at x/L = 3.9. It is seen that in the distance 3.9 = ^- = 5.0, 
this calculation introduces larger errors than those accumulated by the calculation that 
used simple-wave extrapolation all along y/L = 0.75. It is not surprising that the cal- 
culation with the outer boundary at y/L = 1.3 introduces larger errors since, as seen in 
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Figure 22.- Numerical tests on the use of simple-wave extrapolation at 
the outer boundary for a flat-plate flow. 

figure 10, this outer boundary crosses the leading -edge shock -transition region. The 

oscillations in the pressure distribution through the shock in figure 22 are typical of 

those produced by second -order finite -difference schemes. It is clear from the results 

that simple -wave extrapolation cannot be used successfully when the outer boundary is 

placed so that it crosses the shock-transition region. 

Similar tests were conducted with the compression-corner calculations for 
Moo =3.0 and a = 10.0° to determine what differences result because of different posi- 
tions of the outer boundary. Figure 23 shows the four computational boxes for which the 
calculations were made. The outer boundaries of the first two calculations were placed 
between the wall and the leading -edge shock, whereas the outer boundary of the third 
calculation was placed just outside the leading -edge shock at the upstream boundary. The 
outer boundary in the fourth calculation was placed further outside the leading-edge shock, 
and free-stream conditions were imposed up to the approximate position of the intersec- 
tion of the leading -edge shock. Downstream of that point, simple-wave extrapolation was 
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Figure 23.- Computational boxes used to test the simple-wave 
extrapolation for = 3-0 and a = 10.0°. 


used along the outer boundary. Any errors introduced by this extrapolation could not 
propagate upstream since no right-running Mach line intersected the sonic line in the 
boundary layer upstream of the downstream boundary. Since the boundary conditions are 
the most accurate in the fourth calculation, these results should be the most accurate. 
Box I was used for the calculations which were discussed previously. In the first three 
calculations the same grid spacing was used, but in the fourth calculation the following 
variable grid spacing was used in the y -direction in order to reduce the computer 
storage: 

^ = 0.00321 (O ^ i 0.0868) 

x c \ X C / 

expands at 10 percent fo.0868 =^~ = 0.199) > (39) 

x c \ x c / [ 

^ = 0.0136 ( 0.199 § JL < 0.59) 

x c \ x c /_, 

Figure 24 presents the wall pressure distributions which resulted from the four 
calculations. It is seen that the differences shown here are consistent with those found 
in the tests made with the flat -plate calculations. The pressure distribution produced by 
box in is clearly in error since the overall pressure rise is less than the inviscid value. 
The underestimation of the pressure also occurred in the flat -plate tests when simple - 
wave extrapolation was used in the shock -transition region. Without the fourth calcula- 
tion to serve as a reference, it is not as simple to deduce which of calculations one or 
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Figure 2k . - Comparison of wall pressure distributions for three positions 
of the outer boundary of the computational box. 

two is the more correct; hence, this difference is indicative of the uncertainty introduced 
into the calculations by the simple -wave extrapolation at the outer boundary. 

Moo = 4.0 Compression-Corner Calculations 

In this section, calculations are discussed that were made for the Moo = 4.0 flow 
over a 10°. adiabatic compression corner for which Lewis, Kubota, and Lees (ref. 11) 
experimentally obtained wall pressure distributions. The Reynolds number based on 
the free -stream conditions and the distance from the leading edge to the corner was 
Roo,x c = 6.8 x 10*. Considerable care was taken in the experiment to insure that the flow 
was two dimensional and that transition to turbulence occurred downstream of 
reattachment. 

The computational box which was used in the calculations extends from x/x c = 0.38 
to x/x c = 2.11 with the outer boundary placed at a constant distance of 0.109x c from the 
wall. In the x-direction 121 grid points were used and a grid spacing of Ax/x c = 0.01442 
was obtained; in the y-direction 36 grid points were used and a grid spacing of 
Ay/x c = 0.002885 was obtained. An attempt was made to reduce the total number of 
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grid points by increasing Ax by 33 percent; however, the resulting calculation was 
unstable since large oscillations developed in the corner region. 

The upstream boundary conditions were established by a series of calculations from 
the leading edge similar to those discussed previously for the Moo = 3.0 and Moo = 6.06 
flow fields. Details of these calculations are discussed by Carter (ref. 18). The wall 
temperature was held constant at T aw as determined by equation (38) . 

Comparison with experiment .- Comparison of the numerical wall pressure distribu- 
tions with those measured by Lewis et al. are shown in figure 25. The experimental data 
is both from a flat -plate model which spanned the wind tunnel (no side plates) and from 
an axisymmetric model (stovepipe) which consisted of a cylinder -flare configuration. 
These two sets of experimental data agree very well except downstream of the reattach- 
ment region where the axisymmetric pressure distribution begins to drop toward the 
asymptotic cone pressure. 

The agreement between the numerical calculation and the experiment is good 
although downstream of the corner it is seen that the theoretical values are higher. The 
theoretical solution overshoots the inviscid solution; this result is expected based on 


Numerical solution 


Experiment (ref. 11) 



Figure 25. - Comparison of numerical and experimental wall pressure distributions. 
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previous discussion. Calculations were made for two other computational boxes to show 
that the errors introduced by the simple -wave extrapolation are of comparable magnitude 
to those found in the Moo =3.0 calculations previously discussed. These calculations 
also verified that the upstream influence of the corner was contained in the computational 
box. For further details of these calculations, see Carter (ref. 18). 

Figure 26 shows a flow -field map of the streamline and pressure contour pattern 
in the immediate vicinity of the corner. Comparison of figure 26 with the flow -field map 
for the Moo = 3.0 case shown in figure 20 shows that the outer boundary of the computa- 
tional box was placed farther from the wall in the Moo = 4.0 case. The pressure con- 
tours near this outer boundary indicate that no shocks have formed interior to the com- 
putational field. However, in figure 26, the contours are converging and, judging from 
the relatively large distance between the 1.6 and 1.7 pressure contours, it appears that 
distinct separation and reattachment shocks are beginning to form. As was seen in the 
Moo = 3.0 case, the mass flow in the separation bubble is very small compared with that 
in the free-shear layer, although their thicknesses are comparable. 
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Figure 26.- Computed flow field over a 10° compression corner for = 4.0. 


51 



Profiles of flow properties .- In figure 27, profiles of the two velocity components, 
temperature, and density are shown at three stations upstream and downstream of the cor- 
ner. In each case, the profiles are normal to the surface with the velocity components 
parallel and perpendicular to the ramp denoted as u' and v’, respectively. These pro- 
files show the general characteristics of the compression-corner flow field. Figure 27(a) 
shows the large deceleration of the flow near the wall which results in flow separation at 
x/x c = 0.827 as evidenced by the vertical slope of this profile at the wall. At the corner, 
it is seen that the maximum reverse -flow velocity is about 2 percent of the free -stream 
velocity. Also seen in figure 27(a) is the acceleration which occurs in the boundary layer 
after the corner as the flow returns to a normal state of weak interaction at the new Mach 
number. The traverse velocity profiles in figure 27(b) are indicative of the turning of the 
flow. The decrease in v in the outer part of the profile is due to the fact that the outer 
streamlines have not been turned as much as the streamlines closer to the surface. 

Figure 27(c) shows that the temperature in the separation bubble differs only slightly 
from the wall temperature which is expected in view of the low velocities and small nor- 
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(a) u and u 1 component of velocity. 


Figure 27 .- Profiles of flow properties upstream and downstream of the corner. 
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(b) v and. v' component of velocity. 


Figure 27*- Continued. 

mal pressure gradient in this region. The use of a constant wall temperature for the 
approximate adiabatic wall condition is verified by the near -vertical slope of the tem- 
perature profiles in figure 27(c). The wall heat transfer for each profile is less than 
1 percent of the maximum heat transfer in the boundary layer. 

The profiles shown in figure 27 verify, qualitatively at least, the theoretical anal- 
ysis of Stewartson and Williams (ref. 9) which was briefly discussed in the "Introduction." 
A key feature of their analysis is the existence of an inner boundary layer near the sur- 
face where the velocity change is substantial but density changes are small. The inner 
layer is adjacent to the main boundary layer where the u profile is similar in shape to 
that of the unperturbed boundary layer upstream of the ramp. Stewartson and Williams 
(ref. 9) found the thickness of the inner layer to be o(e^) = O^R^Xq) w ^ere Roo,x 0 is 
the Reynolds number based on free-stream conditions and the distance from the leading 
edge to the start of the interaction. In constructing a universal form of the variables for 
the inner region, an amplification of the effective value of e , through a scaling dependent 
on Mach number, results in the thickness of the inner region being of the order 
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Using equation (40) together with the present conditions and x Q /x c = 0.582 results in 
an estimate of y/x c = 0.011 for the inner region thickness. In figure 27(d), the density 
profile at x/x c = 0.827, which is at separation and is within the region of applicability of 
the Stewartson and Williams analysis, shows that this estimate is the correct order of 
magnitude for the lateral extent of the constant -density inner region. Since the pressure 
changes are small, the temperature changes are also small in this region, as seen in fig- 
ure 27(c). Thus, the numerical results tend to verify Stewartson and Williams' use of 
the incompressible boundary -layer equations in this inner region. 

The streamwise velocity perturbation in the main boundary layer was estimated by 
Stewartson and Williams to be of the order 
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(c) Temper ature . 
Figure 2J.~ Continued. 
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which for the present conditions results in Au/u^ ~ 0.27. The amplification of the veloc- 
ity scaling due to the Mach number dependence was small and therefore was deleted. 
Figure 27(a) shows this estimate to be the correct order of magnitude for the difference 
between the profile at separation (x/xc = 0.827) and the profile at the upstream start of 
the interaction (x/xc = 0.582). Similarly, the normal velocity at the outer edge of the 
inner region is of the order 

-3/8 


JL * e 3 = r‘ 

u oo °°> X 0 

whereas in the main boundary layer the normal velocity increases to 
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For R^Xq = 2.93 x 104, equations (42) and (43) result in nondimensional normal veloc- 
ities of the order of 0.02 and 0.07 at the outer edges of inner region and the main boundary 
layer, respectively. This estimate agrees well with the normal -velocity profile at sepa- 
ration shown in figure 27(b). 


|_ = . 582 
x c 



£- = 1.216 

x c 

= 1.432 

x c 

£-=2.h 

x c 



(d) Density. 


Figure 2J.~ Concluded. 
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Moo = 6 .06 Compression-Corner Calculations 

Calculations were made for the Moo = 6.06 flow over a 10.25° adiabatic compres- 
sion corner for which Lewis, Kubota, and Lees (ref. 11) also obtained experimental mea- 
surements for several Reynolds numbers. For Roo,x c = 1.5 X 10^, Lewis et al. judged 
the flow to be laminar and established the repeatability of the surface pressure distribu- 
tion within several percent. Therefore, this case was chosen for the present calculations. 

The computational box used in the calculations extends from x/x c = 0.5 to 
x/x c = 2.1 and from the wall to y/x c = 0.105. The x- and y-grid spacings were 
Ax/x c = 0.0125 and Ay /x c = 0.003, respectively, and resulted in 4644 grid points. It is 
observed that this grid spacing is comparable to that used in the Moo = 4.0 calculation. 
This spacing was used since the boundary -layer thickness is approximately the same in 
both cases because of the compensating effects of increasing the Mach number and 
Reynolds number. The position of the outer boundary of the computational box was far 
enough from the wall to allow for the boundary -layer growth, but did not cross the sepa- 
ration or reattachment shocks determined experimentally from pitot measurements. 

Simple -wave extrapolation was used along the outer boundary except in the region 

0.5 S — = 0.667 where the conditions were known from the flat -plate flow field deter- 

x c 

mined in computational box V which is shown in figure 11. The flow along this part of 
the outer boundary is not influenced by the compression-corner flow field. Tests were 
not made on the simple -wave extrapolation since it was presumed that the errors would 
be comparable to those found previously. 

Figure 28 shows a comparison of the wall pressure distribution determined from 
the present calculations with that found experimentally. Upstream of the corner, the 
present solution and the experiment predict approximately the same degree of compres- 
sion when the upstream pressure difference is taken into account. Downstream of the 
corner, the agreement becomes better as the effects of the upstream pressure difference 
become smaller. Also shown for comparison is the theoretical result obtained by 
Klineberg and Lees (ref. 42) by using an integral method on the boundary -layer equations. 
Their result agrees well with the experiment at the start of the interaction, but down- 
stream it falls below the experimental distribution; as a result, a longer interaction 
region is indicated. 

In figure 29 a comparison between the present results and experiment is shown for 
Mach number profiles at various stations downstream of the corner. In the experiment, 
the boundary -layer thickness was defined on the basis of a pitot profile as is shown in 
figure 15. In order to be consistent, this definition was also used in plotting the pres- 
ent results in figure 29. The agreement obtained in figure 29 is reasonably good, par- 
ticularly at x/x c = 1.2, which is slightly upstream of reattachment, and at x/x c = 2.0; 
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Figure 28.- Comparison of theoretical and experimental wall pressure distributions. 

at x/x c = 1.6 the differences in the experimental and theoretical profiles are greater. 

In the experiment, correction was made where necessary to the pitot measurements 
because of the effects of low Reynolds numbers. 

In figure 30 the skin -friction distribution obtained in the present investigation is 
compared with that found by Klineberg (ref. 43) by using an integral method. The points 
of separation and reattachment predicted by these two approaches differ by about 5 per- 
cent of the. distance from the leading edge. Their locations were not determined experi- 
mentally. Elsewhere, the agreement is reasonably good up to the reattachment region; 
however, downstream of that region, the two distributions differ considerably. The dif- 
ference is consistent with that in figure 28 in that the integral approach shows the inter- 
action to extend farther downstream than that found in the present investigation. 

Murphy (ref. 7) found similar results in comparing the results of three boundary - 
layer approaches (including that of Klineberg and Lees) to experiment. For example, in 
comparing the wall skin friction with that measured by Hakkinen et al. (ref. 23) for the 
case of an incident shock onto a boundary layer for Moo = 2.0, it was found that the extent 
of the separated region was overpredicted by all three methods and that downstream of 
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Figure 29*- Comparison of numerical and experimental Mach number profiles. 

reattachment, the skin friction was underpredicted. Later calculations were made by 
MacCormack (ref. 21) by use of the Navier-Stokes equations for the same case and good 
agreement was obtained with the measured skin friction. 


Comparison With Oswatitsch's Analysis 

Oswatitsch (ref. 44) used Taylor series expansions of the Navier-Stokes equations 
in the vicinity of a separation point to show that the separation streamline (also referred 
to as the dividing streamline) leaves the surface at an angle given by 


0D = tan" 1 



sp /ax 


(44) 
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Present solution 

Integral method 



Figure 30.- Comparison of wall skin-friction coefficient obtained 
from the Navi er -Stokes equations with that from the boundary- 
layer equations. 

In addition, his analysis showed that the locus of u = 0 is a straight line near the sepa- 
ration point inclined to the surface by the angle 


0 U = tan 


-1 


8 T 

9x 
9p/3x 


(45) 


Hence, for small angles, /3 U = ^ j3jy Figure 31 shows a schematic diagram of the stream 
line pattern near separation with the angles /3 U and /3j> indicated. The analysis at 
reattachment is identical if the arrows in figure 31 are reversed. 


Table I shows a comparison between the angles obtained from the calculation and 
those from equations (44) and (45), by use of the present numerical results to estimate 
the shear stress and pressure gradients. With a few exceptions, the differences are less 
than 10 percent; thereby, the results are verified to be consistent with Oswatitsch's 
analysis. 


Free -Interaction Analysis 

Lewis, Kubota, and Lees (ref. 11) correlated their data for the free -interaction 
region by using the scaling laws of Chapman, Kuehn, and Larson (ref. 10) and Curie 
(ref. 45). Curie's scaling is more general since it includes the effect of a nonadiabatic 
wall, but reduces to that of Chapman et al. in the adiabatic case. In Curie's scaling, the 
pressure is given by 
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Dividing streamline 



Figure } 1 .- Schematic diagram of the streamline pattern near 
separation identifying the Oswatitsch angles. 


TABLE I.- COMPARISON OF CALCULATED ANGLES OF DIVIDING STREAMLINE 
AND u = 0 LINE WITH THOSE ESTIMATED BY OSWATITSCH' S THEORY 
FROM THE PRESENT NUMERICAL RESULTS 


Moo 

% 

0u 

Oswatitsch 

Calculated 

Oswatitsch 

Calculated 

Separation 

3.0 

10.6 

10.1 

7.1 

6.5 

4.0 

9.2 

10.0 

6.2 

6.5 

6.06 

10.3 

1 

9.1 

6.9 

6.5 

Reattachment 

3.0 

9.0 

11.0 

6.0 

7.4 

4.0 

7.6 

7.0 

5.1 

4.2 

6.06 

6.0 

6.5 

4.0 

4.1 


60 




(46) 


P = 





and the length scale by 


X cc 




i /4 1/4 

K 0 ,x 0 


(47) 


where the subscript o refers to the flow conditions at x = x Q . The flow conditions at 
x = x Q differ from those at infinity due to the flat -plate weak interaction; hence, Lewis 
et al. used the Lees and Probstein weak-interaction analysis to relate the flow at x = x 0 
to the free -stream values. As a result, equation (46) becomes 


P = 1.752 


P -P 0 1 


0 


1 + 


3y - 1 M +0 ( M del' 

4 °°Vdx/ 0 \ 00 dx 


2 

Jo 


(48) 


and equation (47) becomes, after inserting the same constant of proportionality used by 
Lewis et al., 


X = 0.311 



5 - 3r 
4 



+ 0 M 



(49) 


The variation of the displacement thickness is found from first-order boundary -layer 
theory which gives 


= \jcZ (l.937 ^ + 0.578 




M“ - 0.207 


(50) 


Comparison of the experimental correlation found by Lewis et al. and the present 
results for Moo = 3.0, 4.0, and 6.06 is shown in figure 32. In converting the numerical 
results to the correlation plot, the definition of x 0 used was that employed by Lewis 
(ref. 41), which is shown graphically in that reference. The agreement between the pres- 
ent results and the correlation is good, although the Moo = 6.06 pressure distributions 
are underpredicted by the Moo = 3.0 and Moo = 4.0 correlated results. This difference 
is probably due to the use of a linear pressure -deflection relationship used in deriving 


equation (46), which neglects terms of 


o(m 0 



A second comparison of the present results was made with the universal pressure 
distribution of Stewartson and Williams (ref. 9) for the free -interaction zone. In their 
analysis, the pressure is given by 
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( 51 ) 



and the length scale by 



(52) 


Figure 33 shows the comparison between the present results plotted in these coordinates 
and the numerical distribution obtained by Stewartson and Williams. The X-coordinate 
has been translated so that the separation point is at the origin. The start of the inter- 
action x Q was defined in this comparison as the point upstream of the corner where the 
pressure was a minimum. 


«,deg o,x Q 
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Figure 32.- Comparison of the present numerical results with 
the free- interact ion pressure correlation of Lewis, Kubota, 
and Lees. 




It is seen that the M m = 3.0 and Moo = 4.0 results are well correlated by these 
coordinates and, as in Curie variables, the Moo = 6.06 results do not agree with this 
correlation. This latter result is expected since the Stewartson and Williams free- 
interaction analysis is also based on a linear pressure -deflection model which becomes 
inaccurate in hypersonic flow. Comparison of the present results with that of Stewartson 
and Williams shows that the pressure level is overpredicted by Stewartson and Williams 
by 20 percent. Also, the present results show a slightly greater slope. Stewartson and 
Williams (ref. 9), in comparing with the experimental data of Chapman et al. (ref. 10), 
found differences of comparable magnitude and sign to those found here. The cause of 
the error was unknown, although they suspected that it was due to the fact that their 
asymptotic theory is precise only as e = Since the relative error is 0(e), 

O 

the Reynolds number would have to be 10 in order for this error to be 0.1. 

The correlation variables of Curie differ only in the length scale from those of 
Stewartson and Williams since it is seen that comparing equations (46) and (51) results in 

P 2 = 2" 1/2 P (53) 

In order to compare the x-scaling given in equations (47) and (52), the present results 
were replotted with the origin at the start of the interaction and by using the same defini- 
tion of x 0 , that being the point of minimum pressure. The same degree of correlation 
was obtained with the Curie scaling as that w T ith the Stewartson and Williams scaling; 
thereby, a preference of these two similarity laws could not be made on the basis of the 
present results. 


CONCLUSIONS 

The Brailovskaya finite -difference scheme has been used to obtain steady -state 
solutions to the Navier -Stokes equations for the laminar, supersonic flow over a flat 
plate and a compression corner. In the compression-corner calculations, flows with and 
without separation were considered as well as the effect of wall suction in reducing flow 
separation. Comparison of the calculated results have been made with experimental data 
and with other theoretical solutions. In addition, numerical tests were performed to 
determine the errors introduced by approximate boundary conditions. The following 
conclusions can be made on the basis of the present investigation: 

1. Good agreement between the present results and those from experiment for the 
compression-corner flow field indicates that the Brailovskaya method can be used to obtain 
accurate solutions to the laminar Navier-Stokes equations at high Reynolds numbers. 

2. The good agreement between the solutions of the Navier-Stokes equations and 
similar solutions to the boundary -layer equations indicates that the damping mechanism, 
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which must be present in all explicit schemes in order to maintain stability, is not of 
such an extent for the Brailovskaya scheme that it alters the real viscosity effects. 

3. Caution must be exercised in using approximate techniques (generally some type 
of extrapolation is used) to supply unknown boundary conditions. For example, it was 
found that the use of simple -wave extrapolation in the vicinity of the leading -edge shock 
resulted in erroneous expansion waves being reflected back into the computational region. 

4. The effect of using tpo coarse a grid and neglecting the slip and temperature 
jump boundary conditions in the leading -edge region of a flat plate results only in a 
locally invalid solution. Downstream of this region in the weak-interaction region, good 
agreement was obtained with the weak-interaction theories of Lees and Probstein, and 
Kubota and Ko, except for the skin-friction distribution in the comparison with the analy- 
sis of Kubota and Ko. 

5. Comparison of the present results at a free -stream Mach number of 6.06 with 
those obtained by Klineberg by an integral method of the boundary -layer equations shows 
that the integral method predicts a larger interaction region, particularly downstream of 
the corner. 

6. Comparisons of the calculated angles between the dividing streamline and the 
wall, at separation and reattachment, and those found from equations derived by 
Oswatitsch, agreed, in general, to within 10 percent. 

7. The scaling laws derived by Stewartson and Williams correlated the wall pres- 
sure distributions found in the present calculations up to the corner. The universal 
pressure distribution found numerically by Stewartson and Williams overpredicts the 
pressure level of this correlation by about 20 percent. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., June 1, 1972. 
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APPENDIX A 


STABILITY ANALYSIS OF THE FINITE -DIFFERENCE EQUATIONS 


A von Neumann stability analysis is presented here of the Brailovskaya finite - 
difference scheme applied to the Navier-Stokes equations. This finite-difference scheme 
is explicit and will be conditionally stable provided that it contains a damping mechanism 
which prevents the growth of short wavelength solutions (on the order of a few grid spac- 
ings) to the difference equations. Since the finite -difference scheme is incapable of 
resolving such phenomena, it is necessary for the growth of such solutions to be retarded. 
The von Neumann analysis consists of examining different wavelength solutions to the 
difference equations to establish whether a finite -difference scheme amplifies or damps 
short wavelength solutions. 

This analysis is approximate since only the linearized difference equations are 
used; however, the resulting stability criterion should be viewed as a necessary condition 
since the nonlinearities in the equations of motion would presumably enhance rather than 
retard the amplification of small disturbances. Furthermore, the von Neumann analysis 
does not take into account instabilities which may arise due to the imposition of certain 
boundary conditions. Clearly, the sufficiency of the resulting von Neumann stability 
criterion can only be determined by numerical experimentation. 

As discussed by Richtmyer and Morton (ref. 46), the von Neumann condition for the 
nonconservation form of the governing equations is the same as that for the conservation 
form. Since the analysis is simpler for the nonconservation form, the following non- 
dimensional system of approximate equations, which may be derived from equations (9) 
to (12), were considered: 


-S + A 9 U + b 8 U = e /Aj + AA 

* 9x 9y \9x 2 9y2 J 


where 


(Vl 


u = 


(Al) 
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APPENDIX A — Continued 


A = 


u 

— 

yp 

0 


u 

0 


0 

0 

u 


0 

r - 1 


7 

0 

u 



0 

0 


E = 


0 

4 m/p 
3 Roo 

0 


0 

0 

4 m/p 
3 


0 


0 

0 

0 

yiL 

P 

Np j»Rqo 


and c denotes the speed of sound given by 


c = \J(y - 1)T 


(A2) 


In this system of equations, the viscosity coefficient p 


has been assumed to be constant 
In addition, the mixed deriva- 


in view of the linearization of the von Neumann analysis 
2 2 

fives - v - and - 9 - — ~ have been deleted from the x- and y -momentum equations, 
3x 3y 9x 9y 

respectively, and the dissipation function $ given by 



(A3) 


has been deleted from the energy equation. These last two approximations are based on 
an analysis by Kentzer (ref. 47) which showed that neither mixed derivatives nor terms 
quadratic in first-order derivatives enter into the linear stability consideration. 
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APPENDIX A — Continued 


Application of the Brailovskaya scheme (16) to equation (Al) results in the following 
difference equation where the two steps have been combined by assuming that the coeffi- 
cient matrices A, B, and E are constant: 

°r* 1 ■ °s!k - - "S'-!.*) + m - 2u s.k +u “-2,k) 


AB At (jjii Ijn u n -jj-n \ 

1^4,1 1rii - u ;i1 i - u. i i,.t + k-lj 


, **** [TT aa _ TT 11 _ TT 11 

2 Ax Ay \ 


AE At^ / TT n 


2 Ax 


r( u "+2,k - 2u SVi,k + 2u "-i,k - u ] n -2 >k ) 


AE At^ / TT n 


2 Ax Ay^ 


( U j+1 ,k+l " 2U j+l,k + U j+1 ,k-l " U jbl,k+1 + 2U M,k _ H^-ljk-l) 


B At 
2 Ay 




BE At z / TT n 


2 Ay 3 


• 2 D SW + 2 u SV-i • u ^) 


BE At z / TT n 


2 Ay Ax 


sK + i,k +1 - 2U ",k + l + D f-l,k + X - ^l.k-1 + 2D "k-l - 


E At 
Ax 3 


u" + i,k - + u ”-i,k) ■ < k + u ^-i) 


(A4) 


Substitution of one Fourier component of the solution 
U(x,y,t) = U 0 (t)expjl (k^x + k 2 yjj 

into equation (A4) gives the amplification matrix, which is defined as U?t?= GUl\ 

]>k j,k 


where 


G = 1 - (C 2 + 2D) - iC(I - 2D) 

C _ At (— sln - + B sin £ \ 

V Ax Ay / 


(A5) 

(A6) 

(A7) 
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APPENDIX A - Continued 


D = E At 



cos a 
Ax2 


+ 


1 - cos 
Ay 2 ) 


(A8) 


a - k^ Ax 


0 = k 2 Ay 



(A9) 


and I is the identity matrix. The coefficients kj and k2 denote Fourier components 
of the solution to the difference equations. 

In the present application of the Brailovskaya scheme, calculations are made outside 
of the viscous regions and hence stability must be maintained independent of the magni- 
tude of D. By setting D = 0, equation (A6) becomes 

G = I - C 2 - iC (A10) 


If j u is an eigenvalue of C, the corresponding eigenvalue of G is 

g = 1 - ju 2 - iju (All) 


The von Neumann necessary condition for stability requires that the magnitude of the 
eigenvalues of the amplification matrix be less than one, that is, 

| g | 2 = 1 - M 2 (l - M 2 ) = 1 (A12) 

which will be satisfied provided /i 2 ^ 1. Hence, the maximum eigenvalue of the 
matrix C determines the maximum allowable time increment At. Richtmyer (ref. 27) 
introduced a simple technique by which the eigenvalues of C may be found. Consider 
an axis inclined at an angle 6 with respect to the x-axis where 6 is given by 


cos 6 = 


sin a 
Ax 



(A13) 



(A 14) 


and the velocity component u' along this axis is 
u' = u cos 6 + v sin 6 


(A15) 
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APPENDIX A - Continued 


By use of equations (A13) to (A15), matrix C may be written as 




r* 

u’ 

p cos 6 

p sin 6 




c 2 cos 6 

u’ 

n 

y - 

//sin ol\ 

J /sin /3\ 2 

yp 

U 

y 

|(V Ax ) 

V Ay j 

c 2 sin 6 

0 

u’ 

y - 



yp 



y 



0 

c 2 cos 6 

c 2 sin 8 



cos e\ 
sin 6 


u' 


(A16) 


whose eigenvalues are easily found to be 




* - ♦ (***) L 


lu* + c. 


By inserting equation (A15) for u’ , the stability criterion becomes 

— 2 

1 O 

u sin a v sin 8 


At 


a + ± c 

Ax Ay 


/ sin aV /sin BY 
\ Ax ) \ Ay ) 


Z 1 


which results in the following maximum allowable value for At: 

1 


At S. 


1 \ 2 . / 1 


M + liL + ci/i—i + 

Ax Ay V\ Ax / 


(A17) 


(A18) 


(A19) 


An approximate viscous stability criterion may be derived by neglecting the con- 
vection terms in the Navier-Stokes equations. This approximation is reasonable since, 
as shown by Carter (ref. 18) for Burger’s equation, the viscous stability criterion was 
found to be insensitive to the magnitude of the inviscid terms. Setting C equal to zero 
in equation (A6) results in 

G = I - 2D (A20) 

If g is an eigenvalue of G, then there results 

g = 1 - 2 cp (A21) 
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APPENDIX A - Concluded 


where <p are the eigenvalues of D and are 


r 




<P 


= Atf— - cos a 1 - cos J 
\ Ax 2 Ay 2 / | 


4 m/p 

3 Rco 

4 m/p 
3 Rco 


N p r R °o 
v. 1 ^ 


the maximum of which is the fourth value since for gases 


(A22) 


N 


Pr 


> The von Neumann 

o 


necessary condition for stability requires g SI which results in cp S 1 or 


At S 


O^NpyR^L 



(A23) 


The calculations were made with the minimum At of that given by equation (A19) 
or (A23). In most cases, equation (A19), which is the CFL (Courant -Friedrichs -Lewy) 
condition, gave the smallest At since the Reynolds number was usually large. In those 
cases in which the viscous stability criterion was used, a stable calculation resulted 
thereby suggesting that the approximation made in finding the viscous stability criterion 
is valid. 
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APPENDIX B 


FINITE -DIFFERENCE EQUATIONS ALONG COORDINATE SYSTEM INTERFACE 

Figure 34 shows a typical computational module centered at point 0 along the inter- 
face between a rectangular and skewed coordinate system. The grid spacing is variable 



in the x- and y-directions. Taylor series expansions about the central point to the 
neighboring points results in the following second-order difference expressions: 


d(£ 
9x i 


o Ax 1 AX 2 (AX! + Ax 2 ) 


Ax? (cpr, - Ax., tan a ^ - Ax? tan a JL_S£_ 

2 \ 7 1 ay o 1 ax ay 


l( Ax l tan af 


- Ax^g + (bx\ - Ax|)? 0 


+ 0(A)" 


(Bl) 


where 


d 2 j£_ 


dx dy 
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2 8y 
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9 1 ay 


8 


H - * 4 )^ 
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Ax^ Axg (AXj + AXg^j 


0(42) 


g3 

(b 0 - a Q W 2 - (a Q + b 0 )<p 1 + 2a 0 <p 7 - (b^ - a^) — | 

ay 

ay 9 
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(b Q - a 0 )Ay 1 + (a Q + b Q )Ay 2 


+ 0 


(a 3 ) 
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where 


where 
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APPENDIX B — Continued 


The dependent variable (p represents any of the flow variables p, u, v, or T. In 

are approximated by the usual 


1 

and In. 

V *1 

S’ 3y2 


9y 


three -point difference equations given by equations (22) and (23), respectively. The term 
needs only to be approximated to first order, 


8 3 ^ 

9y 3 


9 


8 3 <p 

9y 3 

To evaluate 

9y 3 I 


_ a 3 <p 

9 8 y 3 


+ 0(A) 


to first order, four points are required, which causes a minor dilemma 


o 


as to which four points to use. The scheme used here is a weighted combination of differ- 
ence expressions obtained from two different staggered sets of points, 

8 3 <p 

8y 3 

where /3 is a weighting factor which varies linearly from zero at the wall to unity at the 
outer boundary of the computational box, and 


= (1 - ft A- 
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(Equations continued on next page) 
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